
PRO MAKE_PHYSICAL_PLOTS_PLOT, x, y, XR=xr, XLOG=xlog, YR=yr, YLOG=ylog, $
                              NBOX=nbox, XBOT=xbot, HIST=hist
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE
  
  ;; Parse Keywords
  xbot = KEYWORD_SET(xbot)
  xlog = KEYWORD_SET(xlog)
  ylog = KEYWORD_SET(ylog)
  
  ;; Make binsize s.t. there are NBOX boxes
  bin1 = xlog ? (alog10(xr[1])-alog10(xr[0]))/nbox : float(xr[1]-xr[0])/nbox
  bin2 = ylog ? (alog10(yr[1])-alog10(yr[0]))/nbox : float(yr[1]-yr[0])/nbox
  print,'Binsize: ',bin1,bin2
  
  ;; Check status of LOG flags and create local variables
  xh  = xlog ? alog10(x)  : x
  xrh = xlog ? alog10(xr) : xr
  yh  = ylog ? alog10(y)  : y
  yrh = ylog ? alog10(yr) : yr
  
  ;; print,'X RANGE: ',xrh
  ;; print,'X Data:  ',m4_stat(xh)
  ;; print,'Y RANGE: ',yrh
  ;; print,'Y Data:  ',m4_stat(yh)
  ;; print,''
  
  ;; Make 2D histogram, and plot
  hist = HIST_2D( xh , yh, $
                  min1 = xrh[0], max1 = xrh[1], bin1 = bin1, $
                  min2 = yrh[0], max2 = yrh[1], bin2 = bin2 )
  
  ;; Load color table, and plot image
  cgLoadct,13,/silent
  plotimage,hist,range=set_plot_range(hist),/noerase,xst=4,yst=4
  
  ;; Set axes so we can see them!
  cgAxis, xaxis=0, /xst, xtickformat='blank_axis', color='wt1', $
          xlog=xlog, xr=xr, /save
  IF ~xbot THEN cgAxis, xaxis=1, /xst, xtickformat='blank_axis', color='wt1'
  cgAxis, yaxis=0, /yst, ytickformat='blank_axis', color='wt1', $
          ylog=ylog, yr=yr, /save
  cgAxis, yaxis=1, /yst, ytickformat='blank_axis', color='wt1'
  
  RETURN
END


PRO MAKE_PHYSICAL_PLOTS, CONFFILE=cfile, NREAL=nreal, NOPLOT=noplot
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE
  
  COMMON OMNI_CONFIG, conf, mw, local, dpdfs, ancil, fmt, conffile
  
  ;; Parse keywords
  noplot = KEYWORD_SET(noplot)
  nreal    = ~n_elements(nreal) ? 1000L : long(nreal)
  
  ;; Read in survey-info, galactic-params, & dpdf-params config files
  conf     = omni_load_conf(cfile)
  mw       = omni_read_conffile('./conffiles/galactic_params.conf')
  dpdfs    = omni_read_conffile('./conffiles/dpdf_params.conf')
  mconf    = omni_read_conffile('./conffiles/mass_deriv.conf')
  d        = dindgen(dpdfs.nbins)*dpdfs.binsize + dpdfs.binstart
  sampletd = n_elements(sampletd) ? sampletd GE 1 : mconf.sampletd
  
  ;; Some physics for a little later...
  c     = 299792458.d                 ; Speed of Light     [m/s]
  h     = 6.626068d-34                ; Planck constant    [m^2 kg s^-1]
  kb    = 1.3806503d-23               ; Boltzmann constant [m^2 kg s^-2 K^-1]
  Jy    = 1.d-26                      ; Jansky             [W m^-2 Hz^-1]
  hnuk  = h*conf.nu/kb                ; h*nu/kb            [K]
  msun  = 1.9891d30                   ; Solar Mass         [kg]
  pc    = 3.08567758d16               ; parsec             [m]
  mu    = 2.8                         ; Avg molecular mass
  mh    = 1.673723599d-24             ; Hydrogen Mass      [g] 
  kappa = conf.kappa/10.d/mw.gas2dust ; dust opactiy       [m^2 kg^-1]
  
  ;; Fiducial mass surface density, modulo flux density (Jy) and size (")
  S0   = c * c * Jy * (exp(hnuk/mw.td) - 1.d) * 206265.d * 206265.d / $
         (2.d * kappa * h * conf.nu*conf.nu*conf.nu*!dpi) * 1.d-1 ; [g cm^-2]
  
  print,'SURFACE DENSITY: ',S0
  
  ;; Create mass functions for use here (sampletd = [0,1])
  message,'Creating SAMPLETD=0 mass function',/inf
  OMNI_MASS_FUNCTION, CONFFILE=cfile, /MEGA, NREAL=nreal, SAMPLETD=0, $
                      /MFNSAVE, MFNFN=fn0, /SILENT
  message,'Creating SAMPLETD=1 mass function',/inf
  OMNI_MASS_FUNCTION, CONFFILE=cfile, /MEGA, NREAL=nreal, SAMPLETD=1, $
                      /MFNSAVE, MFNFN=fn1, /SILENT
  fn = [fn0,fn1]
  tval = ['20K','LogNormal']
  
  restore,fn0
  arsz = size(objs.dist,/DIM)
  nele = n_elements(objs.dist)
  
  flux = dblarr([arsz,2])
  dist = dblarr([arsz,2])
  arad = dblarr([arsz,2])
  tdus = dblarr([arsz,2])
  
  ;; Read in mass function files and compute physical properties
  FOR i=0,1 DO BEGIN
     
     ;; Restore the IDL save file
     restore,fn[i]
     
     ;; Place data into standardized arrays
     flux[*,*,i] = cmreplicate(s.(iflux),n_elements(objs)) ; Flux density [Jy]
     dist[*,*,i] = objs.dist / 1.d3                        ; Distance [kpc]
     arad[*,*,i] = cmreplicate(s.rad,n_elements(objs)) ; Angular rad [arcsec]
     tdus[*,*,i] = temp.objs                           ; Temperature [K]
     
  ENDFOR                        ; End of SAMPLETD loop
  
  
  ;; Compute derived quantities
  radi = arad / 206265.d * (dist *1.d3)                ; Physical radius [pc]
  tcor = (exp(hnuk/tdus)-1.d) / (exp(hnuk/mw.td)-1.d)  ; Temp correction
  mass = m0 * flux * dist * dist * tcor * conf.fluxcor ; M_sun
  dens = mass * msun/ (!dpi*radi*radi*radi*pc*pc*pc) / 1.d3 ; g/cc 
  nden = dens / (mu*mh)                                     ; cm^-3
  sden = S0 * flux / (arad * arad) * tcor                   ; g cm^-2
  cden = sden / (mu*mh)                                     ; cm^-2
  
  ;; Make byte array of elements that have distance in proper range
  usei = dist GE mconf.dmin AND dist LE mconf.dmax AND $
         radi GE mconf.rmin AND radi LE mconf.rmax AND $
         nden GE mconf.nmin AND nden LE mconf.nmax
  fini = finite(arad)
  
  print,'ARAD: '
  print,m4_stat(arad[where(usei)])
  print,m4_stat(radi[where(usei)])
  
  print,''
  print,'Fraction of sources meeting usei criteria:',$
        float(total(usei)) / float(n_elements(usei))
  neari = where(dist LT mconf.dmin AND fini, nnear)
  fari  = where(dist GT mconf.dmax AND fini, nfar)
  print,'Fraction too close: ',float(nnear)/(total(fini))
  print,'Fraction too far: ',float(nfar)/(total(fini))
  
  disti = where(dist GE mconf.dmin AND dist LE mconf.dmax AND fini, ndist)
  radii = where(radi GE mconf.rmin AND radi LE mconf.rmax AND fini, nradi)
  ndeni = where(nden GE mconf.nmin AND nden LE mconf.nmax AND fini, nnden)
  print,ndist,nradi,nnden,nnear,nfar,n_elements(usei)
  
  openw,lun,'./masses/nearfar_limits.txt',/get_lun,/append
  printf,lun,nnear,nfar
  close,lun
  free_lun,lun
  
  IF noplot THEN goto,endofcode
  ;;=====================================================================
  ;; It's plottin' time!
  xsize = 9
  
  
  ;;-----------------------------------------------
  ;; First plot: Mass versus heliocentric distance.
  myps,'./masses/plots/physical_m_d.eps',xsize=xsize
  sun = cgSymbol('sun')
  
  multiplot_xm,[2,1],mpcharsize=1.0
  yr = [3d0,1d4]
  xlog = 0
  ylog = 1
  xtit = 'd'+sun+'  [kpc]'
  
  FOR jj=0,1 DO BEGIN
     
     ;; Extract values for this plot 
     dind = where(usei[*,*,jj]) + nele*jj
     d = dist[dind]
     m = mass[dind]
     
     ;; Setup plot window...
     ytit = ~jj ? 'M  [M'+sun+']' : ''
     ytf  = ~jj ? 'exponent10' : ''
     xr   = [0.d + jj*1d-4, 20.d]
     cgPlot, d, m, /nodata, charsize=1.0, $
             xr = xr, /xst, xlog = xlog, xtit = xtit, $
             yr = yr, /yst, ylog = ylog, ytit = ytit, ytickformat=ytf
     
     ;; Choose number of boxes in each direction, and make plot
     nbox = 100.d               ; 141.4213562d
     MAKE_PHYSICAL_PLOTS_PLOT, d, m, XR=xr, XLOG=xlog, YR=yr, YLOG=ylog, $
                               NBOX=nbox
     
     ;; Add various lines to the plot to help with interpretation 
     vline,mconf.dmin,color='cyan',XLOG=xlog,YLOG=ylog,thick=3
     vline,mconf.dmax,color='cyan',XLOG=xlog,YLOG=ylog,thick=3
     
     ;; Median minimum complete flux density
     dp = findgen(nbox*10 + 1)/(nbox*10) * (xr[1] - xr[0]) + xr[0]
     mp = m0 * 0.28 * dp * dp
     cgOplot,dp,mp,color='bisque',thick=3,linestyle=3
     
     ;; Label lines
     cgText,14,500,'S!d1.1!n = 0.28 Jy',ali=0.5,chars=0.9,col='bisque',orie=22
     
     al_legend,/bottom,/right,/clear,['T!dd!n = '+tval[jj]]
     
     multiplot,/doxaxis,doyaxis=0
  ENDFOR
  myps,/done,/mp
  
  
  
  
  
  ;;-----------------------------------------------
  ;; Second plot: Mass versus Physical Radius.
  myps,'./masses/plots/physical_m_r.eps',xsize=xsize
  sun = cgSymbol('sun')
  
  multiplot_xm,[2,1],mpcharsize=1.0
  yr = [3d0,1d4]
  xlog = 1
  ylog = 1
  xtit = 'R  [pc]'
  
  FOR jj=0,1 DO BEGIN
     
     ;; Extract values for this plot 
     dind = where(usei[*,*,jj]) + nele*jj
     r = radi[dind]
     m = mass[dind]
     
     ;; Setup plot window...
     ytit = ~jj ? 'M  [M'+sun+']' : ''
     ytf  = ~jj ? 'exponent10' : ''
     xr   = [4d-2 + jj*1d-4, 1.d1]
     cgPlot, r, m, /nodata, charsize=1.0, $
             xr = xr, /xst, xlog = xlog, xtit = xtit,$
             yr = yr, /yst, ylog = ylog, ytit = ytit, ytickformat=ytf
     
     ;; Choose number of boxes in each direction, and make plot
     nbox = 100.d               ; 141.4213562d
     MAKE_PHYSICAL_PLOTS_PLOT, r, m, XR=xr, XLOG=xlog, YR=yr, YLOG=ylog, $
                               NBOX=nbox
     
     ;; Add various lines to the plot to help with interpretation 
     vline,0.125,XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     vline,1.25, XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     
     ;; Constant Density lines for Dunham(2011) definitions [750,1d4]
     rp = 10.d^( findgen(nbox*10 + 1)/(nbox*10) * $
                 (alog10(xr[1]) - alog10(xr[0])) + alog10(xr[0]) )
     mp = 750.d * (mu*mh) * 1.d3 / msun * (pc*pc*pc) * (rp*rp*rp) * !dpi
     cgOplot,rp,mp,color='bisque',thick=3,linestyle=3
     mp = 1.d4 * (mu*mh) * 1.d3 / msun * (pc*pc*pc) * (rp*rp*rp) * !dpi
     cgOplot,rp,mp,color='bisque',thick=3,linestyle=3
     
     ;; Label regions for core/clump/cloud
     cgText,.07,6,'Cores',align=0.5,color='wt1',charsize=0.9
     cgText,4.5,500,'Clouds',align=0.5,color='wt1',charsize=0.9
     cgText,0.30,15,'Clumps',align=0.5,color='wt1',charsize=0.9
     
     res = linfit(alog10(r),alog10(m))
     print,'LINFIT RES: ',res
     xpl  = 10.^(!x.crange)
     ypl  = 10.^(res[0])*xpl^(res[1])
     ypl2 = 228.*xpl^(2.36)
     cgOplot,xpl,ypl,color='brown',thick=5,linestyle=4
     cgOplot,xpl,ypl2,color='chartreuse',thick=5,linestyle=4
     al_legend,/top,/left,color=['brown','chartreuse'],thick=5,linestyle=4,$
               linsize=0.6,box=0,textcolor='wt1',$
               [string(10.^res[0],res[1],$
                       format="('y = ',F0.1,' x!u',F0.2,'!n')"),$
                'y = 228 x!u2.36!n']
     
     al_legend,/bottom,/right,/clear,['T!dd!n = '+tval[jj]]
     
     multiplot,/doxaxis,doyaxis=0
  ENDFOR
  myps,/done,/mp
  
  
  
  
  
  ;;-----------------------------------------------
  ;; Third plot: Mass versus Surface Density
  myps,'./masses/plots/physical_m_sd.eps',xsize=xsize
  sun = cgSymbol('sun')
  
  multiplot_xm,[2,1],mpcharsize=1.0
  yr = [3d0,1d4]
  xlog = 1
  ylog = 1
  xtit = cgSymbol('Sigma')+'  [g cm!u-2!n]'
  
  FOR jj=0,1 DO BEGIN
     
     ;; Extract values for this plot 
     dind = where(usei[*,*,jj]) + nele*jj
     sd = sden[dind]
     m  = mass[dind]
     
     ;; Setup plot window...
     ytit = ~jj ? 'M  [M'+sun+']' : ''
     ytf  = ~jj ? 'exponent10' : ''
     xr   = [1.d-3 + jj*1d-6, 1.d0]
     cgPlot, sd, m, /nodata, charsize=1.0, $
             xr = xr, /xst, xlog = xlog, xtit = xtit, xtickformat='exponent10',$
             yr = yr, /yst, ylog = ylog, ytit = ytit, ytickformat=ytf
     
     ;; Choose number of boxes in each direction, and make plot
     nbox = 100.d               ; 141.4213562d
     MAKE_PHYSICAL_PLOTS_PLOT, sd, m, XR=xr, XLOG=xlog, YR=yr, YLOG=ylog, $
                               NBOX=nbox
     
     ;; Add various lines to the plot to help with interpretation 
     vline,0.0107,XLOG=xlog,YLOG=ylog,color='bisque',thick=3,linestyle=3
     ;; vline,1.d4, XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     
     ;; Constant Radius lines for Dunham(2011) definitions [0.125,1.25]
     np = 10.^(findgen(1001)/1000*3 - 3.) * 1.d1 / msun
     mp = !dpi * (0.125 * pc)^2 * np
     cgOplot,np/1.d1*msun,mp,color='bisque',thick=3,linestyle=2
     mp = !dpi * (1.25 * pc)^2 * np
     cgOplot,np/1.d1*msun,mp,color='bisque',thick=3,linestyle=2
     
     ;; Constant nDen lines for Dunham(2011) definitions [750,1d4]
     mp = !dpi * np*np*np / (750.d * 750.d * (mu*mh*mu*mh) * 1.d6 / msun/msun)
     cgOplot,np/1.d1*msun,mp,color='bisque',thick=3,linestyle=2
     mp = !dpi * np*np*np / (1.d4 * 1.d4 * (mu*mh*mu*mh) * 1.d6 / msun/msun)
     cgOplot,np/1.d1*msun,mp,color='bisque',thick=3,linestyle=2
     
     
     ;; Label regions for core/clump/cloud
     cgText,3d-1,1.2d1,'Cores',align=0.5,color='wt1',charsize=0.9
     cgText,2.5d-3,5.d3,'Clouds',align=0.5,color='wt1',charsize=0.9
     cgText,8.d-3,6,'Clumps',align=0.5,color='wt1',charsize=0.9
     
     ;; Label lines
     cgText,0.3,85,'R = 0.125 pc',chars=0.9,col='bisque',ali=0.5,orie=42
     cgText,2.5d-3,75,'R = 1.25 pc',chars=0.9,col='bisque',ali=0.5,orie=42
     cgText,3.7d-3,10,'n = 750 cm!u-3!n',chars=0.9,col='bisque',ali=0.5,orie=69
     cgText,0.16,1d3,'n = 10!u4!n cm!u-3!n',chars=0.9,col='bisque',$
            ali=0.5,orie=69
     cgText,9.5d-3,3.7d3,'2$\sigma$ rms noise',chars=0.9,col='bisque',$
            ali=0.5,orie=90
     cgText,1.45d-2,3.7d3,'@T!dd!n = 20K',chars=0.9,col='bisque',$
            ali=0.5,orie=90
      
     
     
     al_legend,/bottom,/right,/clear,['T!dd!n = '+tval[jj]]
     
     multiplot,/doxaxis,doyaxis=0
  ENDFOR
  myps,/done,/mp
  
  
  
  
  
  ;;-----------------------------------------------
  ;; Fourth plot: Mass versus Number Density
  myps,'./masses/plots/physical_m_n.eps',xsize=xsize
  sun = cgSymbol('sun')
  
  multiplot_xm,[2,1],mpcharsize=1.0
  yr = [3d0,1d4]
  xlog = 1
  ylog = 1
  xtit = 'n  [cm!u-3!n]'
  
  FOR jj=0,1 DO BEGIN
     
     ;; Extract values for this plot 
     dind = where(usei[*,*,jj]) + nele*jj
     n = nden[dind]
     m = mass[dind]
     
     ;; Setup plot window...
     ytit = ~jj ? 'M  [M'+sun+']' : ''
     ytf  = ~jj ? 'exponent10' : ''
     xr   = [5d1 + jj*1d-1, 4.d5]
     cgPlot, n, m, /nodata, charsize=1.0, $
             xr = xr, /xst, xlog = xlog, xtit = xtit,$
             yr = yr, /yst, ylog = ylog, ytit = ytit, ytickformat=ytf
     
     ;; Choose number of boxes in each direction, and make plot
     nbox = 100.d               ; 141.4213562d
     MAKE_PHYSICAL_PLOTS_PLOT, n, m, XR=xr, XLOG=xlog, YR=yr, YLOG=ylog, $
                               NBOX=nbox
     
     ;; Add various lines to the plot to help with interpretation 
     vline,750,XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     vline,1.d4, XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     
     ;; Constant Radius lines for Dunham(2011) definitions [0.125,1.25]
     np = 10.^(findgen(1001)/1000*5 + 1.)
     mp = !dpi * (0.125 * pc)^3 * np * mu * mh / msun * 1.d3
     cgOplot,np,mp,color='bisque',thick=3,linestyle=2
     mp = !dpi * (1.25 * pc)^3 * np * mu * mh / msun * 1.d3
     cgOplot,np,mp,color='bisque',thick=3,linestyle=2
     
     ;; Constant Surface Density Lines
     mp = !dpi * 1.0^3 / (np*np*mu*mu*mh*mh) * (1.d-3/msun)
     cgOplot,np,mp,color='bisque',thick=3,linestyle=3
     mp = !dpi * 0.1^3 / (np*np*mu*mu*mh*mh) * (1.d-3/msun)
     cgOplot,np,mp,color='bisque',thick=3,linestyle=3
     mp = !dpi * 0.01^3 / (np*np*mu*mu*mh*mh) * (1.d-3/msun)
     cgOplot,np,mp,color='bisque',thick=3,linestyle=3
     
     ;; Label regions for core/clump/cloud
     cgText,1.5d5,10,'Cores',align=0.5,color='wt1',charsize=0.9
     cgText,1.2d2,6.d3,'Clouds',align=0.5,color='wt1',charsize=0.9
     cgText,1.9d3,4,'Clumps',align=0.5,color='wt1',charsize=0.9
     
     ;; Label lines
     cgText,1.5d5,75,'R = 0.125 pc',chars=0.9,col='bisque',ali=0.5,orie=50
     cgText,150,75,'R = 1.25 pc',chars=0.9,col='bisque',ali=0.5,orie=50
     
     cgText,1.7d5,3.5d3,cgSymbol('Sigma')+' = 1.0 g cm!u-2!n',ali=0.5,$
            chars=0.9,col='bisque',orie=-67
     cgText,2d4,2.5d2,cgSymbol('Sigma')+' = 0.1 g cm!u-2!n',ali=0.5,chars=0.9,$
            col='bisque',orie=-67
     cgText,2.5d3,1.7d1,cgSymbol('Sigma')+' = 0.01 g cm!u-2!n',ali=0.5,$
            chars=0.9,col='bisque',orie=-67
     
     al_legend,/bottom,/right,/clear,['T!dd!n = '+tval[jj]]
     
     multiplot,/doxaxis,doyaxis=0
  ENDFOR
  myps,/done,/mp
  
  
  
  
  
  ;;-----------------------------------------------
  ;; Fifth plot: Number density versus heliocentric distance
  myps,'./masses/plots/physical_n_d.eps',xsize=xsize
  sun = cgSymbol('sun')
  
  multiplot_xm,[2,1],mpcharsize=1.0
  yr = [5d1,4d5]
  xlog = 0
  ylog = 1
  xtit = 'd'+sun+'  [kpc]'
  
  FOR jj=0,1 DO BEGIN
     
     ;; Extract values for this plot 
     dind = where(usei[*,*,jj]) + nele*jj
     n = nden[dind]
     d = dist[dind]
     
     ;; Setup plot window...
     ytit = ~jj ? 'n  [cm!u-3!n]' : ''
     ytf  = ~jj ? 'exponent10' : ''
     xr   = [0.d + jj*1d-4, 20.d]
     cgPlot, d, n, /nodata, charsize=1.0, $
             xr = xr, /xst, xlog = xlog, xtit = xtit, $
             yr = yr, /yst, ylog = ylog, ytit = ytit, ytickformat=ytf
     
     ;; Choose number of boxes in each direction, and make plot
     nbox = 100.d               ; 141.4213562d
     MAKE_PHYSICAL_PLOTS_PLOT, d, n, XR=xr, XLOG=xlog, YR=yr, YLOG=ylog, $
                               NBOX=nbox
     
     ;; Add various lines to the plot to help with interpretation 
     vline,mconf.dmin,color='cyan',XLOG=xlog,YLOG=ylog,thick=3
     vline,mconf.dmax,color='cyan',XLOG=xlog,YLOG=ylog,thick=3
     vline,/h,750,XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     vline,/h,1.d4, XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     
     ;; Constant Radius lines for Dunham(2011) definitions [0.125,1.25]
     np = 10.^(findgen(1001)/1000*5 + 1.)
     mp = !dpi * (0.125 * pc)^3 * np * mu * mh / msun * 1.d3
     cgOplot,np,mp,color='bisque',thick=3,linestyle=2
     mp = !dpi * (1.25 * pc)^3 * np * mu * mh / msun * 1.d3
     cgOplot,np,mp,color='bisque',thick=3,linestyle=2
     
     ;; Label regions for core/clump/cloud
     cgText,1.5d5,10,'Cores',align=0.5,color='wt1',charsize=0.9
     cgText,2.d2,4.d3,'Clouds',align=0.5,color='wt1',charsize=0.9
     cgText,3.d3,6,'Clumps',align=0.5,color='wt1',charsize=0.9
     
     al_legend,/bottom,/right,/clear,['T!dd!n = '+tval[jj]]
     
     multiplot,/doxaxis,doyaxis=0
  ENDFOR
  myps,/done,/mp
   
  
  
  
  
  ;;-----------------------------------------------
  ;; Sixth plot: Physical Raduis versus heliocentric distance
  myps,'./masses/plots/physical_r_d.eps',xsize=xsize
  sun = cgSymbol('sun')
  
  multiplot_xm,[2,1],mpcharsize=1.0
  yr = [4d-2,1d1]
  xlog = 0
  ylog = 1
  xtit = 'd'+sun+'  [kpc]'
  
  FOR jj=0,1 DO BEGIN
     
     ;; Extract values for this plot 
     dind = where(usei[*,*,jj]) + nele*jj
     r = radi[dind]
     d = dist[dind]
     
     ;; Setup plot window...
     ytit = ~jj ? 'R  [pc]' : ''
     ytf  = ~jj ? 'exponent10' : ''
     xr   = [0.d + jj*1d-4, 20.d]
     cgPlot, d, r, /nodata, charsize=1.0, $
             xr = xr, /xst, xlog = xlog, xtit = xtit, $
             yr = yr, /yst, ylog = ylog, ytit = ytit, ytickformat=ytf
     
     ;; Choose number of boxes in each direction, and make plot
     nbox = 100.d               ; 141.4213562d
     MAKE_PHYSICAL_PLOTS_PLOT, d, r, XR=xr, XLOG=xlog, YR=yr, YLOG=ylog, $
                               NBOX=nbox
     
     ;; Add various lines to the plot to help with interpretation 
     vline,mconf.dmin,color='cyan',XLOG=xlog,YLOG=ylog,thick=3
     vline,mconf.dmax,color='cyan',XLOG=xlog,YLOG=ylog,thick=3
     vline,/h,0.125,XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     vline,/h,1.25, XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     
     ;; Physical size associated with BGPS spatial transfer function
     dp = findgen(nbox*10 + 1)/(nbox*10) * (xr[1] - xr[0]) + xr[0]
     rp = 14.d / 206265. * dp * 1.d3
     cgOplot,dp,rp,color='bisque',thick=3,linestyle=3
     ;; rp = 33.d / 206265. * dp * 1.d3
     ;; cgOplot,dp,rp,color='bisque',thick=3,linestyle=3
     rp = 150.d / 206265. * dp * 1.d3
     cgOplot,dp,rp,color='bisque',thick=3,linestyle=3
     
     ;; Label lines
     cgText,14,0.75,cgSymbol('theta')+'!dR!n = 14"',ali=0.5,chars=0.9,$
            col='bisque',orie=16
     cgText,5,4,cgSymbol('theta')+"!dR!n = 2.5'",ali=0.5,chars=0.9,$
            col='bisque',orie=34
     
     al_legend,/bottom,/right,/clear,['T!dd!n = '+tval[jj]]
     
     multiplot,/doxaxis,doyaxis=0
  ENDFOR
  myps,/done,/mp
  
  
  
  
  
  ;;-----------------------------------------------
  ;; Seventh plot: Physical Raduis versus number density
  myps,'./masses/plots/physical_n_r.eps',xsize=xsize
  sun = cgSymbol('sun')
  
  multiplot_xm,[2,1],mpcharsize=1.0
  yr = [3d1,3d5]
  xlog = 1
  ylog = 1
  xtit = 'R  [pc]'
  
  FOR jj=0,1 DO BEGIN
     
     ;; Extract values for this plot 
     dind = where(usei[*,*,jj]) + nele*jj
     r = radi[dind]
     n = nden[dind]
     
     print,'Median Values: '
     print,'n: ',10^(m4_stat(alog10(n)))
     print,'R: ',10^(m4_stat(alog10(r)))
     
     ;; Setup plot window...
     ytit = ~jj ? 'n  [cm!u-3!n]' : ''
     ytf  = ~jj ? 'exponent10' : ''
     xr   = [4d-2 + jj*1d-4, 1.d1]
     cgPlot, r, n, /nodata, charsize=1.0, $
             xr = xr, /xst, xlog = xlog, xtit = xtit, $
             yr = yr, /yst, ylog = ylog, ytit = ytit, ytickformat=ytf
     
     ;; Choose number of boxes in each direction, and make plot
     nbox = 100.d               ; 141.4213562d
     MAKE_PHYSICAL_PLOTS_PLOT, r, n, XR=xr, XLOG=xlog, YR=yr, YLOG=ylog, $
                               NBOX=nbox
     
     ;; Add various lines to the plot to help with interpretation 
     vline,0.125,XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     vline,1.25, XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     vline,750.,/h,XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     vline,1.d4,/h,XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     
     ;; Various Surface Density Contours
     rp = findgen(nbox*10 + 1)/(nbox*10) * (xr[1] - xr[0]) + xr[0]
     np = 1.0 / (mu*mh) / (rp * pc * 1.d2)
     cgOplot,rp,np,color='bisque',thick=3,linestyle=3
     np = 0.1 / (mu*mh) / (rp * pc * 1.d2)
     cgOplot,rp,np,color='bisque',thick=3,linestyle=3
     np = 0.01 / (mu*mh) / (rp * pc * 1.d2)
     cgOplot,rp,np,color='bisque',thick=3,linestyle=3
     
     ;; Larson's Relationships
     larcol = 'deep pink'
     np = 3420. * (2.d * rp)^(-1.10d)
     cgOplot,rp,np,color=larcol,thick=5,linestyle=4
     np = 840. * (2.d * rp)^(-0.51d)
     cgOplot,rp,np,color=larcol,thick=5,linestyle=4
     
     cgText,8d-2,1.8d5,'Cores',align=0.5,color='wt1',charsize=0.9
     cgText,2.5d+0,1d2,'Clouds',align=0.5,color='wt1',charsize=0.9
     cgText,2.0d-1,6d3,'Clumps',align=0.5,color='wt1',charsize=0.9
     
     ;; Label Lines
     cgText,4,2.2d4,cgSymbol('Sigma')+' = 1.0 g cm!u-2!n',ali=0.5,chars=0.9,$
            col='bisque',orie=-32
     cgText,3,3d3,cgSymbol('Sigma')+' = 0.1 g cm!u-2!n',ali=0.5,chars=0.9,$
            col='bisque',orie=-32
     cgText,0.25,1.6d3,cgSymbol('Sigma')+' = 0.01 g cm!u-2!n',ali=0.5,$
            chars=0.9,col='bisque',orie=-32
     
     cgText,0.07,3.6e4,'Larson',ali=0.5,chars=0.9,col=larcol,orie=-36
     
     al_legend,/bottom,/right,/clear,['T!dd!n = '+tval[jj]]
     
     multiplot,/doxaxis,doyaxis=0
  ENDFOR
  myps,/done,/mp
  
  
  
  
  
  ;;-----------------------------------------------
  ;; Eighth plot: Radius versus Surface Density
  myps,'./masses/plots/physical_r_sd.eps',xsize=xsize
  sun = cgSymbol('sun')
  
  multiplot_xm,[2,1],mpcharsize=1.0
  yr = [4d-2,1d1]
  xlog = 1
  ylog = 1
  xtit = cgSymbol('Sigma')+'  [g cm!u-2!n]'
  
  FOR jj=0,1 DO BEGIN
     
     ;; Extract values for this plot 
     dind = where(usei[*,*,jj]) + nele*jj
     sd   = sden[dind]
     m    = mass[dind]
     sd2  = sd * 10. / msun * pc * pc
     
     ;; Setup plot window...
     ytit = ~jj ? 'R  [pc]' : ''
     ytf  = ~jj ? 'exponent10' : ''
     xr   = [1.d-3 + jj*1d-6, 1.d0]
     cgPlot, sd, r, /nodata, charsize=1.0, $
             xr = xr, /xst, xlog = xlog, xtit = xtit, xtickformat='exponent10',$
             yr = yr, /yst, ylog = ylog, ytit = ytit, ytickformat=ytf
     
     ;; Choose number of boxes in each direction, and make plot
     nbox = 100.d               ; 141.4213562d
     MAKE_PHYSICAL_PLOTS_PLOT, sd, r, XR=xr, XLOG=xlog, YR=yr, YLOG=ylog, $
                               NBOX=nbox, /XBOT
     cgAxis, xaxis=1, /xst, color='black',charsize=1.0,$
             xr=10.^(!x.crange)*10./msun*pc*pc,xtickformat='exponent10'
     cgAxis, xaxis=1, /xst, color='wt1',charsize=1.0,$
             xr=10.^(!x.crange)*10./msun*pc*pc,xtickformat='blank_axis'
     ;; Add various lines to the plot to help with interpretation 
     vline,0.0107,XLOG=xlog,YLOG=ylog,color='bisque',thick=3,linestyle=3
     ;; vline,1.d4, XLOG=xlog,YLOG=ylog,color='bisque',thick=3
     
     ;; ;; Constant Radius lines for Dunham(2011) definitions [0.125,1.25]
     ;; np = 10.^(findgen(1001)/1000*3 - 3.) * 1.d1 / msun
     ;; mp = !dpi * (0.125 * pc)^2 * np
     ;; cgOplot,np/1.d1*msun,mp,color='bisque',thick=3,linestyle=2
     ;; mp = !dpi * (1.25 * pc)^2 * np
     ;; cgOplot,np/1.d1*msun,mp,color='bisque',thick=3,linestyle=2
     
     ;; ;; Constant nDen lines for Dunham(2011) definitions [750,1d4]
     ;; mp = !dpi * np*np*np / (750.d * 750.d * (mu*mh*mu*mh) * 1.d6 / msun/msun)
     ;; cgOplot,np/1.d1*msun,mp,color='bisque',thick=3,linestyle=2
     ;; mp = !dpi * np*np*np / (1.d4 * 1.d4 * (mu*mh*mu*mh) * 1.d6 / msun/msun)
     ;; cgOplot,np/1.d1*msun,mp,color='bisque',thick=3,linestyle=2
     
     
     ;; ;; Label regions for core/clump/cloud
     ;; cgText,3d-1,1.2d1,'Cores',align=0.5,color='wt1',charsize=0.9
     ;; cgText,2.5d-3,5.d3,'Clouds',align=0.5,color='wt1',charsize=0.9
     ;; cgText,8.d-3,6,'Clumps',align=0.5,color='wt1',charsize=0.9
     
     ;; ;; Label lines
     ;; cgText,0.3,85,'R = 0.125 pc',chars=0.9,col='bisque',ali=0.5,orie=42
     ;; cgText,2.5d-3,75,'R = 1.25 pc',chars=0.9,col='bisque',ali=0.5,orie=42
     ;; cgText,3.7d-3,10,'n = 750 cm!u-3!n',chars=0.9,col='bisque',ali=0.5,orie=69
     ;; cgText,0.16,1d3,'n = 10!u4!n cm!u-3!n',chars=0.9,col='bisque',$
     ;;        ali=0.5,orie=69
     cgText,9.5d-3,0.09,'2$\sigma$ rms noise',chars=0.9,col='bisque',$
            ali=0.5,orie=90
     cgText,1.45d-2,0.09,'@T!dd!n = 20K',chars=0.9,col='bisque',$
            ali=0.5,orie=90
      
     
     
     al_legend,/bottom,/right,/clear,['T!dd!n = '+tval[jj]]
     
     multiplot,/doxaxis,doyaxis=0
  ENDFOR
  myps,/done,/mp
  
  
  
  
  

  
  
  
  
  
  
  endofcode:
  RETURN
END
